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Monolayer graphene is an example of materials with multi-valley electronic structure. In such 
materials, the valley index is being considered as an information carrier. Consequently, relaxation 
mechanisms leading to loss of valley information are of interest. Here, we calculate the rate of valley 
relaxation induced by charged impurities in graphene. A special model of graphene is applied, where 
the Pz orbitals are two-dimensional Gaussian functions, with a spatial extension characterised by an 
effective Bohr radius OeB- We obtain the valley relaxation rate by solving the Boltzmann equation, 
for the case of noninteracting electrons, as well as for the case when the impurity potential is screened 
due to electron-electron interaction. For the latter case, we take into account local-field effects and 
evaluate the dielectric matrix in the random phase approximation. Our main findings: (i) The valley 
relaxation rate is proportional to the electronic density of states at the Fermi energy, (ii) Charged 
impurities located in the close vicinity of the graphene plane, at distance d < 0.3 A, are much more 
efficient in inducing valley relaxation than those farther away, the effect of the latter being suppressed 
exponentially with increasing graphene-impurity distance d. (iii) Both in the absence and in the 
presence of electron-electron interaction, the valley relaxation rate shows pronounced dependence on 
the effective Bohr radius Obb- The trends are different in the two cases: in the absence (presence) of 
screening, the valley relaxation rate decreases (increases) for increasing effective Bohr radius. This 
last result highlights that a quantitative calculation of the valley relaxation rate should incorporate 
electron-electron interactions as well as an accurate knowledge of the electronic wave functions on 
the atomic length scale. 

PACS numbers: 71.45.Gm, 72.10.Fk, 72.80.Vp, 76.20.+q 


I. INTRODUCTION 

Certain crystalline solids, such as monolayer and 
bilayer graphene, carbon nanotubes, transition-metal 
dichalcogenides, silicon, and diamond, possess multi¬ 
valley electronic structure. Recently, ways to con¬ 
trol and measure the valley degree of freedom (or 
valley index, for short) in these materials have been 

proposecPll^ and tested experimentallyvalley in¬ 
dex is also being actively considered as a carrier of quan¬ 
tum information.ESinHlll 

The valley index of an electron is linked to its crys¬ 
tal momentum. The crystal momentum is changed upon 
scattering, and consequently, the electron can be moved 
between different valleys by scattering processes {inter¬ 
valley scattering) in an uncontrolled, random fashion. 
Therefore, such scattering processes lead to the loss of 
information encoded in the valley index. 

In this work, we theoretically study how (classical) val¬ 
ley information, encoded in an ensemble of electrons, is 
lost due to scattering processes. We focus on a specific 
multi-valley material, monolayer graphene ,^211111 because 
of its relatively simple band structure and widespread 
experimental availability. In particular, we consider in¬ 
travalley and intervalley scattering of graphene’s elec¬ 
trons off nearby charged impurities (Coulomb scatter¬ 
ing), and calculate the corresponding valley relaxation 
time, that is, the time scale characterizing the loss of 
valley information. 

Coulomb scattering is a well-studied mechanism 


as a determinant of the electrical conductivity of 
graphene,^2lH3l] q^j. knowledge, its role in inter¬ 

valley scattering has not been studied in detail. In fact, 
the Coulomb potential of a charged impurity is expected 
to be less efficient in inducing intervalley scattering than 
in inducing intravalley scattering, since the Fourier spec¬ 
trum of the Coulomb potential is peaked around small 
wave numbers. Here, we set out to go beyond that quali¬ 
tative argument by quantifying the efficiency of Coulomb 
scattering for valley relaxation. 

At present, the understanding of valley relaxation 
processes is rather limited; the various mechanisms, 
such as electron-phonon, electron-impurity and electron- 
electron scattering, and their material-specific details, 
are yet to be systematically investigated. Note, how¬ 
ever, that recent studies have started to elucidate vari¬ 
ous aspects of valley relaxation and decoherence in two- 
dimensional (2D) transition-metal dichalco genidePH^ 
and graphenepl^ as w ell as carbon-basecP^^ and sili¬ 
con quantum dots.l^^^^Sl 

To describe valley relaxation in graphene due to 
charged impurities, we use a special model of graphene’s 
electrons, in which the orbitals are described by 
2D Gaussian functions, with a spatial extension char¬ 
acterised by an effective Bohr radius Oes ■ We obtain the 
valley relaxation rate Tv by solving the corresponding 
Boltzmann equation, for the case of noninteracting elec¬ 
trons, as well as for the case when the impurity potential 
is screened due to electron-electron interaction. For the 
latter case, we calculate the screened impurity potential 
by taking into account local-field effects and evaluating 
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the dielectric matrix in the random phase approximation 

(RPA)S3B3| 

Our main findings are as follows, (i) The valley re¬ 
laxation rate is proportional to the electronic density 
of states at the Fermi energy, (ii) Charged impurities 
located in the close vicinity of the graphene plane, at 
distance d < 0.3 A, are much more efficient in inducing 
valley relaxation than those farther away, the effect of 
the latter being suppressed exponentially with increasing 
graphene-impurity distance d. (iii) Both in the absence 
and in the presence of electron-electron interaction, the 
valley relaxation rate shows pronounced dependence on 
the effective Bohr radius OeB • Remarkably, the trends are 
different in the two cases: in the absence (presence) of 
screening, the valley relaxation rate decreases (increases) 
for increasing effective Bohr radius. This last result high¬ 
lights that a quantitative calculation of the valley relax¬ 
ation rate should incorporate electron-electron interac¬ 
tions as well as an accurate knowledge of the electronic 
wave functions on the atomic length scale. 

It should be emphasized that intervalley scattering has 
consequences beyond inducing valley relaxation: it has 
its fingerprints on the magnetoconductivity as well as on 
inelastic light scattering, i.e., the Raman spectrum. In 
graphene, the quantum correction to the conductiv ity is 
influenced by elastic intervalley scattering processesPSESl 
for weak (strong) intervalley scattering, the correction 
to the conductivity is positive (negative), corresponding 
to weak antilocalisation (weak localisation). In experi¬ 
ments, a negative correction can be observed which is at¬ 
tributed to a significant intervalley scattering rate.^^sHUl 
In Raman spectra, the D peak intensity increases with 
increasing intervalley scattering.!^ Furthermore, the in¬ 
tervalley scattering rate can be monitored in real space 
via spatially resolved Raman spectroscopy,!^ revealing 
that the sample boundary can be a strong source of in¬ 
tervalley scattering. 


II. PRELIMINARIES 

We define the valley polarization as the imbalance 
of the electronic populations in the two valleys {riK and 
n^/), i.e., riv = uk — n-K'- Our aim is to describe the 
dynamics of the valley-polarized initial state shown in 
Fig. under the influence of impurity scattering. It 
is expected that impurity scattering transfers electrons 
from one valley to the other (intervalley scattering), and 
therefore leads to the decay of the valley polarization 
with time, see Fig. [^. The task is to quantify the time 
evolution of this decay. 

The qualitative nature of the dynamics of the electron 
distribution depends strongly on the relative time scales 
of elastic and inelastic scattering process. In this work, 
we will describe how elastic intervalley scattering caused 
by static impurities contributes to the decay of valley 
polarization; hence, we will disregard inelastic processes. 
Under the assumption that inelastic processes are absent. 



FIG. 1. Evolution of a valley-polarized initial state due to 
scattering processes, (a) Initial state with a finite valley po¬ 
larization riv (b) Non-equilibrium valley-unpolarized state, 
which appears during the relaxation process if the elastic in¬ 
tervalley scattering is faster than the inelastic processes, (c) 
Thermal equilibrium, reached from (a) or (b), due to inelas¬ 
tic scattering processes, (d) Schematic representation of the 
time evolution of valley polarization. The characteristic time 
scale of the decay is the valley relaxation time Tv. 


the valley-polarized initial state shown in Fig. evolves 
to the non-equilibrium, but valley unpolarized state de¬ 
picted in Fig. [^. If inelastic intravalley transitions are 
present, but they are slow compared to the elastic in¬ 
tervalley processes, then they reinforce thermal equilib¬ 
rium (Fig. 0=) after the intermediate state in Fig. is 
reached. However, if inelastic intravalley transitions are 
faster than elastic intervalley processess, then the initial 
state of Fi g. [T^ evolves directly toward thermal equilib¬ 
rium (Fig. IR). Even though we do not incorporate in¬ 
elastic processes in our model below, we expect that our 
treatment provides an accurate description of the valley 
relaxation time in all cases discussed above. 

It is customary to distinguish short-range and long- 
range impurities. The short-range label usually refers 
to crystallographic defects, charge-neutral adatoms, etc. 
Long-range refers to charged scatterers that induce a 
long-range Coulomb potential for the mobile electrons. 
These impurities might be located, e.g., in the substrate 
supporting the graphene sheet, as shown in Fig. For 
the present paper, we consider a model where the relax¬ 
ation of valley polarization is due to charged impurities 
that are randomly positioned in a plane at a given dis¬ 
tance d from the graphene sheet (see Fig. [^). (Gen¬ 
eralization of our methods to other impurity types and 
spatial distributions is probably straightforward.) 

Figure shows a single charged impurity located at 
a distance d from the graphene sheet. Assuming this is a 
negatively charged impurity with charge — e, it creates 
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FIG. 2. Charged impurities as a source of valley relax¬ 
ation. (a) Random spatial arrangement of charged impurities 
(red spheres) in the substrate (gray) supporting the graphene 
sheet, (b) Schematic representation of the impurity-induced 
Coulomb potential in the plane of the graphene sheet. 


a repulsive potential energy landscape Vi{r;d) for the 
delocalized electrons in graphene; here r = (x, y) is the 
position vector in the graphene plane. If the impurity is 
located on the z axis of the reference frame, then the 2D 
Fourier transform of the potential V) reads 

yi(q;d) = ^e-«^ (1) 

where we use Cq = e^/47reo. Here, e is the magnitude of 
the electron charge and cq is the vacuum permittivity. 


III. MODEL 

A. Boltzmann equation and the valley relaxation 
time 

In this section, we calculate the valley relaxation rate 
due to elastic electron-impurity scattering in the frame¬ 
work of Boltzmann theory. We find that the valley po¬ 
larization riv decays exponentially with time, and obtain 
a simple relation [Eq. (0] between the corresponding 
rate, i.e., the valley relaxation rate and the momentum- 
dependent intervalley scattering rates. 

For clarity, we assume that the valence band is filled, 
and the distribution function characterizing the occupa¬ 
tion of bulk conduction-band states is denoted by /fc, 
i.e., it is assumed to be independent of position. Then 
the Boltzmann equation reads: 

( 2 ) 

h' 

where W^k' is the impurity-induced transition rate from 
state k to state k'. Here we assumed detailed balance 


Wkk' = bFfc'fci which is reasonable as the scattering rates 
will be evaluated using Fermi’s Golden Rule: 

StT \ 

^^kk' ~y~ I ^ck^ck' I ^ (^cfc ^ck ') ■ (11) 

Here, Vck.ck' = {ck\V\ck') is the matrix element of the 
impurity potential V{r;d) = ~ with 

the Bloch-type energy eigenstates |cfc) of the conduc¬ 
tion band (see below). The impurity potential Vlr\d) 
is the electric potential energy created by the impurity 
ensemble with impurities. The overline denotes dis¬ 
order average, with the impurity positions ri assumed 
to be independently and homogeneously distributed. Sck 
denotes the energy of the conduction-band state at mo¬ 
mentum k. 

The solution of the Boltzmann equation will be facili¬ 
tated by the fact that in graphene, the low-energy excita¬ 
tions have an approximately conical (linear and i sotropic) 
dispersion relation around the Dirac points.l^SlIIlHere, we 
consider a sample which has a Fermi energy e-p that is 
small compared to the energy width of the tt band, the 
latter being approximately 16.8 eV. This has the follow¬ 
ing implications. First, the dispersion is approximately 
conical (see Fig. at energy ep, and the two Fermi lines 
are approximately circles with radii Kp = ep/hvp, where 
Up ~ 9 X 10® m/s is the Fermi velocity. Second, the 
crystal momenta k of the electronic states participat¬ 
ing in the valley relaxation process are close to either 
K or K' , hence we can uniquely relabel their momen¬ 
tum ask = K + K or k = K' -b k, respectively, with 
K = |k| < K. 

Let us briefly and qualitatively discuss the energy 
range where the conical (linear and isotropic) approx¬ 
imation of the electronic dispersion is valid. We have 
checked that the error of the conical approximation [see 


Eq. (18)] with respect to the tight-binding dispersion 
[see Eq. (16l] is at most 10% for wave vectors fulfilling 


K < Kc = 0.316/acc = 2.23 x 10® m i.e., within this 
wave-vector range the relations 0.9 < en,K+K < 

fulfilled. The Fermi energy corresponding 
to Kc is hvpKc ~ 1.33 eV. Therefore, it is reasonable to 
use the conical approximation of the electronic dispersion 
as long as 


ep < 1.33 eV. 


(4) 


The Boltzmann-type description of carrier dynamics 
in graphene in the presence of Coulomb scatterers is ex¬ 
pected to be valid if the carrier density exceeds a thresh¬ 
old set by the graphene-impurity distance d and the im¬ 
purity density n, = Ni/A, where A is the sample area. 
The theoretical and experimental grounds of this expec¬ 
tation are summarized in, e.g., Sec. HI. A of Ref. [SI] 
and in Ref. ini There are strong indications that the 
Boltzmann-type description is invalid for very low carrier 
densities, when the Fermi energy is in the close vicin¬ 
ity of the Dirac points; see, e. g., Refs. 1351 and IMl and 
Sec. IV. of Ref. 
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The zero-temperature thermal equilibrium state of 
the electrons in this sample is described by the zero- 
temperature Fermi-Dirac distribution /o(£cfc) = ©(^f ~ 
Ecfe). The initial state we consider is depicted in Fig. &■ 
The distribution function in the initial state (Fig. can 
be formulated using valley-dependent Fermi energies Efk 
and £fk'- 


Jk+k (t = 0) = Q{£fk — £c,k+k.), (5a) 

Ik'+k (t = 0) = Q{£fK' — £c,K'+k)- (5b) 

In line with Fig. we assume sfk' < £f < £fk', due to 
particle conservation, only two of the three Fermi ener¬ 
gies are independent. Furthermore, we consider an initial 
state where the occupation difference in the two valleys 
is small, in the sense that all the scattering rates of the 
states participating in the dynamics can be approximated 
by scattering rates at the Fermi energy ep. 

Our aim here is to describe the relaxation dynamics of 
the valley polarization density n^{t), which is related to 
the distribution function via 


with 


Fv = 




( 11 ) 


provided that the sum on the right hand side of Eq. 0 
is independent of the direction of n. This latter condi¬ 
tion is fulfilled in the case we consider, and this becomes 
apparent when evaluating the integral in Eq. (26) below. 

Combining Eqs. §), 0, and ( [l0| ), we find that the 
valley polarization density also shows an exponential de¬ 
cay: 


nv(t) = nv(0)e 


-Fvt 


( 12 ) 


Therefore, we call Fv the valley relaxation rate and Tv = 
Fv ^ the valley relaxation time. Note that according to 
Eq. @, the valley relaxation rate is twice as large as 

the intervalley scattering rate \ 


B. Dispersion, wave functions and scattering rates 


’^v(0 = ^0[/t^+i^(i)-/*:'+ic(^)], (6) 

where A is the sample area. This relation suggests that 
in order to obtain n^{t), it is not necessary to solve the 
original Boltzmann equation ([^. Instead, formulating 
and solving a time-evolution equation for the distribution 
difference 


= Ik+k - Ik'+k. (7) 

might be sufficient. This can be done if the conditions 


Wk+u.k+k.' = lEj 


K'+n.K'+n' = (8a) 

WK+n,K’+n' = Wk'+^.K+k.' = (8b) 


are fulfilled. Eqs. (Sal and (8b) describe the intravalley 
and intervalley transition rates, respectively. From now 
on, we rely on these conditions, and in Appendix we 
argue that they are indeed approximately fulfilled under 
the small-Fermi-energy condition [Eq. 

A straightforward calculation using Eqs. (|^, (^3 and 
(8b) shows that the time-evolution equation for reads 

(v) 


dr. 


dt 


= (/I)’-A') 

k' 


( 9 ) 


Importantly, the initial distribution difference (0), 
dehned via Eqs. 0 and is approximately isotropic in 
n in due to our small-Fermi-energy condition [Eq. @]- 
Therefore the Boltzmann equation (|^ is solved by t 
time-evolving distribution function 


ne 


(i) = /^^ (0)< 


.-Fvi 


( 10 ) 


In the previous section, we established the relation be¬ 
tween the valley relaxation time Tv and the intervalley 
scattering rates ^ ■ The latter can be obtained from 

Fermi’s Golden Rule in Eq. (|^, if the Bloch-type elec¬ 
tronic wave functions '4’ckii') and the dispersion relation 
£ck of the conduction band are known. Here, we outline 
the model we use to evaluate these quantities. 

We use the standard tight-binding or LCAO (linear 
combination of atomic orbitals) model of the tt electrons 
of graphene, with a special feature that the atomic Pz 
orbitals are represented by normalized 2D Gaussian-like 
wave functions, 

0(r) = —=1-e '‘“eB, (13) 

V ^TTOeB 

where OeB is the characteristic length scale of the or¬ 
bitals, and the label eB corresponds to ‘effective Bohr 
radius’. We use OeB as a parameter and assume that 
it is approximately an order of magnitude smaller than 
the carbon-carbon distance. Following Ref. [HU we will 
disregard the overlap of atomic wave functions centered 
on different atoms. This is a reasonable approximation 
as long as OeB < 0.3 acC) as the latter inequality guaran¬ 
tees that the nearest-neighbor overlap integral is less than 
0.25. (For a free-standing carbon atom, we estimate that 
the spatial extension of a three-dimensional 2^^ atomic 
orbital in the xy plane is approximately 0.28 accj see 
Appendix [d|) 

With the above simplifications, the electronic Bloch 
wave functions are expressed as 

Ink) = ^ e*''^ [onk \(t>R) + bnk |^|)) , (14) 

where is the number of the unit cells in the sample, 
(^I0 r) = — R) and {r\cj)^) = 4>{r — R — r), A and 











5 


B are the two sublattices, ii is a lattice vector, and 
b„k are the sublattice amplitudes that can be obtained 
from the tight-binding modeP^^ as 

te) - (t) 

with f(fe) = 1-1- e“‘^“i -I- e“'^“^. Here, the band index 
n e (- 1 - 1 ,— 1 ) = (c,v) refers to the conduction (- 1-1 or c) 
and valence (—1 or v) bands. Note that the vector r, 
and the primitive lattice vectors ai and a 2 , are defined 
in Appendix]^ and the direct lattice, the reference frame 
and the reciprocal lattice are shown in Fig.j^of Appendix 
[Al The dispersion relation reads 

£nk = nsk = n 7 o|f(fc)|, (16) 

where 70 ~ 2.8 eV is the nearest neighbour hopping ma¬ 
trix element. 

In the vicinity of the two Dirac points, f{k) can be lin¬ 
earized, yielding the approximate sublattice amplitudes 


I ^nK+n 

1 

\^nK+K./ 

^nK' -\-n 
\ nK' -\-n / 


1 / ne-'A 
V2[ 1 ) 

1 ) 


(17a) 

(17b) 


where ip = Z {k, x) is the polar angle of k, and the lin¬ 
earized dispersion relation 


= ^nK'+n = (18) 


where up = 37oacc/2^ « 9 x 10® m/s is the Fermi veloc¬ 
ity in graphene. 

Using the 2D Gaussian-like atomic wave function (13) 
in our tight-binding model is a simplification. In princi¬ 
ple, one could use more realistic models, e.g., the linear 
combination of three-dimensional hydrogen-type atomic 
orbitals, or Bloch wave functions obtained from a numeri¬ 
cal density-functional calculation. As we show below, our 
choice (131 has the advantage that it yields simple ana¬ 
lytical expressions for the calculated quantities, includ¬ 
ing the inverse of the dielectric matrix at wave vector K. 
Furthermore, although our simplified approach might not 
yield quantitatively accurate results, it is expected to re¬ 
veal the qualitative role of the atomic wave function in 
the intervalley scattering processes, and to be used as a 
benchmark for future numerical approaches. 


IV. RESULTS 

A. Unscreened impurities 

Here, we consider the case when the screening of the 
impurity Coulomb potential due to electron-electron in¬ 
teraction is disregarded. Using the scattering rate in 


Eq. (§ and the isotropic linear spectrum around the 
Dirac points in Eq. (18), we obtain a formula from 
Eq. for the valley relaxation rate 


y. _ 2eFA 


- 

/ c\ \ ^C,JC-f 
JO 
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k’ —K—K^ 


(19) 


where p' = Z. {k', x) is the polar angle of k'. 

The scattering matrix element can be obtained using 
the electronic wave function in Eq. (14). By neglecting 


the contributions of integrals involving atomic wave func¬ 
tions located at different sites, we find 

^ 'y ^ ^c,K-\-k,c,K' +K,' T G) 

G 

X P {K + Ak + G)V* {K + Ak + G; d), 

( 20 ) 

where we use that K' — K is equivalent with K on the re¬ 
ciprocal lattice. (Note that K = 47r/3-\/3acc ~ 1.7 A 
with acc ~ 1.42 A being the carbon-carbon distance 
in the graphene lattice.) Furthermore Am = k' — k, 
P {q) = / d^r e““*'’|(() (r)l^ = e“2“=B9^ jg the Fourier- 
transformed probability density of the atomic orbital 
(‘form factor’), V {q;d) = J d^re~''^''V {r;d) is the 
Fourier-transformed impurity potential and 

^nk^n'h' (^) ^nk^n'k' 4“ ^nk^n'k'^ ^ (21) 

is a factor from sublattice amplitudes (‘structure factor’). 
We have also used that P{q) is real-valued, which is im¬ 
plied by the cylindrical symmetry of the atomic wave 
function, (j)[r) = (f) (r). 

Because of our small-Fermi-energy assumption 
[Eq. (H] , the condition k,k' ^ K holds and therefore 
Ak K. That implies that Eq. (20) can be well 
approximated by taking the limit Am —>■ 0: 




c,K~{-k,c,K 




j(lin) 


, {K + G) 


G 


X P{K+ G)V* (K+ G;d). (22) 

Here, we kept m and m' in the lower index of the struc¬ 
ture factor; using the approximated sublattice ampli¬ 
tudes Eqs. (17a) and (17b) we obtain 




a-i{K+G)T _ pi(<^+v') 


-.(23) 


The next step is to perform the disorder average in 


Eq. (19). Note that V{q;d) = Vi{q] and 

the assumption of homogeneously and independently po¬ 
sitioned impurities implies 


V* (K + G)V(K + G') = NiV^ (K + G; d) Sgg' ■ 


(24) 


Using Eq. ( |24[ ), we find 

rEEEEEf = § E (k + G] 


X P‘^{K + G) {K + G-,d). (25) 
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Using Eq. (231 to evaluate the integral over the polar 
angle ip' of Kq yields 


^271- 


dLp' 

27r 


c(lin) 

‘^c,K+k.,c,K'+k,‘ 


, {K + G) 


1 

2 ’ 


(26) 


which implies that the valley relaxation rate reads 
ni£F 


Ev = 




^P'^{K + G) [K + G;d). (27) 


G 


Here, ni = N-JA is the sheet density of the impurities. 
Equation (27) is one of the three key results of this work. 


Let us now discuss the parameter dependence of the 
valley relaxation rate Ev. First, we note that its depen¬ 
dence on TT-i is linear as expected. Second, the valley 
relaxation rate Ev is proportional to the Fermi energy 
Ef- This linear relation is a simple consequence of the 
fact that the number of final states in the K' valley, that 
can be reached by scattering from the K valley, is pro¬ 
portional to the density of states D{e-p) = 
per spin per valley, which, in graphene, is proportional 
to e-p itself. The property that the valley relaxation rate 
is proportional to the density of states at the Fermi en¬ 
ergy is expected to be true for other multi-valley mate¬ 
rials as well. For example, in 2D monolayer transition- 
metal dicalcogenides, where the low-energy dispersion re¬ 
lation in the conduction band is parabolic, we expect a 
Fermi-energy-independent valley relaxation rate (for un¬ 
screened impurities). 


Another characteristic feature of the result (271 is the 
G sum, where squared Fourier components of the im¬ 
purity potential are summed up with the weight func¬ 
tion . The appearance of the Fourier components 
Vi {K + G) is not surprising, since in our model, the bulk 
wave functions |cfe) are not plane waves but Bloch-type 
wave functions, i.e., superpositions of plane waves with 
wave numbers fe - 1 - G. Note also that the number of rel¬ 


evant terms in the G sum of Eq. (27) is controlled by 
the effective Bohr radius OeB- If the effective Bohr ra¬ 
dius is increased, the momentum-space weight function 


P\q) = 


a-OeB9 


becomes narrower, and the number of 


Fourier components that contribute significantly to the 
valley relaxation rate decreases. Consequently, the val¬ 
ley relaxation time grows with increasing effective Bohr 
radius. 

This trend is seen in Fig. where the valley relax¬ 
ation time corresponding to the unscreened result (27), 


for graphene-impurity distance d = 0 , is shown as the 
blue line. The result was obtained by numerically com¬ 
puting the sum in Eq. ( |27[ ), which converges due to the 
Gaussian decay of P^(fcb T 


of Fig. 


Eq. (27), i.e., 


The unit on the vertical axis 
^ is defined as the G = 0 term of the sum in 


v,o(c^) — 




n -1 


2 t^2 

“ 2Kd 


dTr^CQUi^F 


(28) 


where Eq. 0 was used, furthermore, Ev,o(d) = Tvo(d). 
For example, using the realistic parameter set rii = 



FIG. 3. Valley relaxation time as a function of the effective 
Bohr radius. The charged impurities are assumed to be lo¬ 
cated in the graphene plane (d = 0). The cases of unscreened 
impurities (blue) and screened impurities (red) are shown. 
The green curve corresponds to a screened calculation where 
the off-diagonal matrix elements of the invrese dielectric ma¬ 
trix are disregarded (‘diagonal screening’). The unit of the 
vertical axis is defined in Eq. (28l. The unit of the horizontal 
axis is the carbon-carbon distance acc- 


10^^ cm“^, Ef = 0.1 eV and d = 0, we find rv,o ~ 10 ps. 
Note that the corresponding transport lifetime for the 
same parameter set is Ttr ~ 5 fs. 

For a given d, the value of rv,o(d) can be regarded 
as an order-of-magnitude estimate of the valley relax¬ 
ation time. Equation (28) reveals that this estimate 
grows exponentially with graphene-impurity distance d, 
Tv^o oc where the characteristic length scale of the 
growth is £ = l/ 2 Ar « 0.3 A. This approximately ex¬ 
ponential dependence of the valley relaxation time, and 
the above-estimated characteristic length scale, are illus¬ 
trated in Fig. 1^, where Ev is shown as the function 
of the graphene-impurity distance d. One implication 
of this approximately exponential behaviour is as fol¬ 
lows. If the graphene layer is lying directly on a sub¬ 
strate, then charged impurities in the latter might be at 
a few-angstrom distance from the graphene plane, and 
can cause a relatively short valley relaxation time. For 
example, if rii = 10 ^^cm“^, ef = 0.1 eV and d = 3 A, 
then rv,o ~ 200 ns. However, if the graphene layer is 
suspended at a finite distance above the substrate, then 
the valley relaxation time improves exponentially with 
the distance; e.g., if d = 3nm, then Tv,o improves with a 
factor of « 2 X 10"^. 

A further thing to note about the G sum of Eq. (27) is 
that the three terms corresponding to G = 0, G = — bi 
and G = 62 are equal (see Fig. for the definitions of bi 
and 62 )) because the corresponding three wave vectors 
K, K — bi and K + b 2 are equal in length, and Vi{q) 
and P{q) are both cylindrically symmetric. In fact, these 
three terms dominate the G sum in the case of large 
graphene-impurity distance, d > acG: as the remain¬ 
ing terms corresponding to longer wave vectors are sup- 



























7 




FIG. 4. Valley relaxation rate as a function of graphene- 
impurity distance (unscreened impurities). (a) Approxi¬ 
mately exponential decay of the valley relaxation rate with 
graphene-impurity distance. Valley relaxation rate rv(d) is 
shown for three different values of the effective Bohr radius. 
The unit of the vertical axis is Fv^o at d = 0 [see Eq. (281]. 
(b) Deviations from exponential decay: same data as in (a), 
but here the ratio of the valley rel axa tion rate and the d- 
dependent rv,o(d) oc [see Eq. (28l] is plotted. The plot 

indicates that for d < acc, the decay of the valley relaxation 
rate with d does not follow oc exactly, and the devi¬ 

ation becomes less significant as the effective Bohr radius is 
increased. 


on this result by taking into account electron-electron in¬ 
teraction, which screens the impurity-induced Coulomb 
potential. In this subsection, we use the RPA approach 
to evaluate the dielectric matrix of graphene at the in¬ 
tervalley wave number, then we use that to 

calculate the valley relaxation time in Sec. |IV C| 

We note that most of the preceding theoretical works 
studying dielec tric screening by the tt electrons of 
graphen^^SEHm apply jellium-type descriptions, where 
the charge density corresponding to the Bloch-type wave 
functions is assumed to be homogeneous, the electronic 
dispersion relation is approximated by the Dirac cones, 
the two valleys are treated independently, their contri¬ 
butions to screening are simply added up, and screening 
is described by a dielectric function e{q) instead of the 
dielectric matrix we use below. This method might be 
appropriate as long as one is interested in the screening of 
the long-wave-length Fourier components of the impurity 
potential. This is the case, e.g., when the conductivity is 
calculated, since that is largely determined by intraval¬ 
ley scattering processes. However, here we consider in¬ 
tervalley scattering, which involves a momentum trans¬ 
fer comparable to the inverse of the atomic length scale. 
Therefore, it is required that the atomic-scale structure 
of the electronic wave functions of the crystal, as well 
as the atomic-scale structure of the electrostatic poten¬ 
tial (a. k. a. local-field effects), are taken into account, 
and that the wave-vector summations are performed for 
the Brillouin zone. We do these following the approach 
of Adlei!^ and Wiser.l^ We note that the dielectric re¬ 
sponse of graphene has been described, with local-field 
effects included, to characterise intervalley plasmons,!^ 
and the macroscopic static dielectric function.l^ 

To characterise how the impurity Coulomb potential 
is screened by the electrons in graphene, we need to ex¬ 
press the total (or screened) potential Vtot that has two 
contributions: the external potential I4xt created by the 
impurities, and the induced potential created by the rear¬ 
ranged electrons. This relation is customarily expressed 
via the inverse dielectric function e~^{r,r') as 

Vfcot(»’)= j d^r (r,r') Dext (r-')- (29) 


pressed due to the relation V^{q; d) oc see Eq. (0. 

As a consequence, Ty{d) oc rv,o(d) holds for d > acc- 
This relation is demonstrated in Fig.[^, where the ratio 
Tv(d)/Tv o(d) is indeed shown to saturate for d > acc- 
Figure also shows that, in the case oi d < acc, the 
proportionality relation Tv oc rv,o breaks down, especially 
for smaller values of the effective Bohr radius CeB- 


B. Screening within the random phase 
approximation 

So far, we have evaluated the valley relaxation time in 
the presence of unscreened impurities. Next, we improve 


Graphene, being crystalline, has discrete translational in¬ 
variance. This implies that an external potential with 
wave vector q induces a total potential which is a super¬ 
position of Fourier components with wave vectors q + G, 
where G is a reciprocal lattice vector. Hence the rela¬ 
tion between the Fourier components of the external and 
induced potential reads 

(q + G) = J2 ^GG' (9) ^ext (q + G') , (30) 

G' 

where q is defined in the first Brillouin zone and G, G' 
are reciprocal vectors. The quantity e~^[q) is called the 
inverse dielectric matrix, its matrix elements are denoted 
by ef^Q,{q), and these matrix elements are related to the 


























Fourier components of the inverse dielectric function via 


^GG' ( 9 ) = e 


-1 


G,-q-G') 


(31) 


We emphasise that the inverse dielectric matrix is a 
matrix-valued function whose domain is the first Bril- 
louin zone. The dielectric matrix e{q) is also a matrix¬ 
valued function on the first Brillouin zone, fulfilling 
[e(q)]“^ = e~^{q) for all q. 

The dielectric matrix is related to the polarizability 
matrix n(< 7 ) as (see Appendix [C|): 


The intervalley polarizability matrix n(lT) can be ex¬ 
pressed by invoking Eqs. (34), (33), ( [^ , and (15), re¬ 
spectively. The result is 

Hgg' (K) = UkP iK + G)P{K + G') [l + e‘(®'-®) 

(35) 

where we defined 

_ ffs / (Snk) ~ / (Sn'k+K) (36) 

4A ^ Enk — Sn'k+K 
nn' k 


sgg' ( 9 ) — ^GG' - Vc (9 + G) IIgg/ (q), 


(32) Note that the quantity f is absent from the result (|35 ), 


where 5gg' is the Kronecker delta, and Vq [q + G) = 
27reg/|q + G| is the 2D Fourier transform of the Coulomb 
potential of the electron-electron interaction. 

In the RPA, the polarizability matrix is approximated 
with that of the noninteracting electron system. The 
latter can be obtained from first-order static perturb ation 
theory, and is expressed via the Adler-Wiser formulcPSEl 


ncG^ (9) = f E 


/ (^nfc) / {_^n'k+q) 


ements in Eq. (33) can be simplified to 


(nfc|e + q) = P{q + G) Snk.n'k+q {q + G), 

(34) 

where P {q + G) and Snk,n'k' (9 + G) are defined above. 

Recall that our goal is to use the inverse dielectric ma¬ 
trix for calculating the valley relaxation rate. To this 
end, we need to know the Fourier components of the to¬ 
tal potential in the vicinity of the wave vectors K + G 
only. Therefore, we need to evaluate the inverse dielec¬ 
tric matrix in the vicinity of K. Assuming that the 
inverse dielectric matrix is a smooth function, we will 
calculate e~^{K) explicitly and use the approximation 
e~^{K -t- Ak) « e~^{K) whenever K + Ak is close to 
K or an equivalent wave vector. Note that by using the 
notation e~^{K) above, we have implicitly defined K to 
be part of the first Brillouin zone. Furthermore, we will 
call e{K), and n(iT) the intervalley dielectric 

matrix, the inverse intervalley dielectric matrix, and the 
intervalley polarizability matrix, respectively. 


even though f appears in Eq. (15) describing the sublat¬ 


tice amplitudes; the reason is that the k sum [in Eq. (33)] 
of the terms containing f vanishes. 

The quantity tljc is evaluated assuming zero temper¬ 
ature T = 0 and charge neutrality ep = 0 implying 
f{e) — 0(—e). I.e., the conduction band is empty and 
the valence band is fully occupied, hence only the in¬ 
terband (n n') terms contribute to the n, n' sum in 
Eq. ([3^: 


(33) 


^nk ^n'k-\-q 

X -|- q) {n'k + q|e'('^^® ('’jnfc) , 

where gs = 2 accounts for the twofold spin degener¬ 
acy. Note that Eq. ( [33| ) is a generalisation of the Lind- 
hard formula.^ The latter expresses the polarizability of 
a homogeneous noninteracting system, whereas the for¬ 
mer generalises that to the case of an inhomogeneous 
system with discrete translational invariance. Further¬ 
more, as we describe the regime of small Fermi ener¬ 
gies, we will calculate and use the polarizability matrix 
of charge-neutral graphene (at zero temperature), which 
corresponds to f{e) = 0(—e) in Eq. (33). 

Neglecting the contributions of integrals involving 
atomic wave functions at different sites, the matrix el- 




£k + Sk+K 


(37) 


where Sk was defined in Eq. (16). To evaluate the sum¬ 


mation over k, we first convert it to a momentum-space 
integral. Despite the ^ 1/k divergence of the inte¬ 
grand, the integral is well defined due to its 2D nature. 
Numerical integration yields Hk ~ —0.143/700^^ « 
-2.53 X 10-2l/eVA^ 

The intervalley dielectric matrix e{K) can be ob¬ 
tained using Eq. ( |3^ and Eq. ( [35| ). To determine the 
screened potential, we need to find the inverse inter¬ 
valley dielectric matrix e~^{K), fulfilling the relation 






-I 


V [K) = 5c 


We find that the 


matrix elements of e ^(iT) are 


^GG' (-^) = ^GG' + 


Vc{K + G)T\gg'{K) 


Sk 


(38) 


where we introduced the dimensionless quantity 


= 1 - E {K^G). (39) 


G 


The analytical formula (38) for the inverse intervalley 


dielectric matrix is the second one of the three key results 
of this work. 

Let us close this subsection by discussing the qualita¬ 
tive features of the results. 

First, focus on the intervalley polarizability matrix 
n(lT) given in Eq. (35). (1) The characteristic scale of 
its matrix elements is given by the quantity !!/<-, which 
is independent of the effective Bohr radius OeB- (2) The 
diagonal elements of n(lT) are real; however, the off- 
diagonal elements are complex in general, because of the 
factor in the square brackets. (3) The dependence of 
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nGG'(-f^) on OeB can be interpreted as follows. First, 
recall that the weight function P{q) is a 2D Gaussian 
function with a characteristic momentum-space width of 
1 /oeB ■ Therefore, Eq. (321 testifies that an external po¬ 
tential of wave vector K + G' can effectively create an 
induced electron density of wave vector K + G, ii and 
only if both wave vectors are below 1 /UeB ■ If either K+G 
OT K + G' are longer than l/oeB, then either P{K + G) 
or P{K + G') suppresses the matrix element nGG'(-I^)- 
This makes sense: On the one hand, if the wave vec¬ 
tor K + G' characterising the external potential is much 
longer then l/ogB, then the corresponding wave length 
is much shorter than UeBj hence the effect of the poten¬ 
tial on the electronic wave functions ‘averages out’ and is 
therefore small indeed. On the other hand, the induced 
electron density is composed of atomic orbitals, hence 
it cannot accommodate spatial variations with smaller 
wave length than OeB- Therefore its Fourier spectrum 
is constrained to the wave vectors shorter than l/ceBj 
explaining the suppression factor P{K + G) in Eq. (32). 


Second, we consider the quantity eji relevant for the 
inverse intervalley dielectric matrix. By approximating 
the sum in Eq. (39) to an integral, we find 


ejt « 1 - IIj^ 


d^fc 


Ibi X br 


(fe) (fe) « 1 + 


1.19 


OeB/aCC’ 
(40) 

where | bi x b 2 1 is the momentum-space area of the Bril- 
louin zone. In Fig. we show ck as a function the effec¬ 
tive Bohr radius, as obtained via numerical evaluation of 
Eq. (391 (solid line) and via the integral approximation in 
Eq. (401 (dashed line); the two results show a reasonable 
qualitative agreement. 



FIG. 5. Effective intervalley dielectric constant as a function 
of the effective Bohr radius. The unit of the horizontal axis 
is the carbon-carbon distance acc- The results of numerical 
summation and analytical approximation are shown. 


C. Screened impurity 


Here, we use our result (381 for the inverse dielectric 


matrix to calculate the valley relaxation rate correspond¬ 
ing to screened charged impurities. I.e., the disordered 
potential V{r-d) appearing in the scattering rate ([^ is 
assumed to take the form [cf. Eq. (29)] 


Ni 


V{r-,d)= J d^r'e ^(r, r') ^ Ei(r'— d). (41) 

i=i 


Repeating the calculation presented in Sec. |IV A| with 
this disorder potential, we find 


Fv = 


niSF 


^ 1_Lpi(G'-G)^ 

E - P{K + G)P{K + G') 


GG'G" 


X fee" (^)]' 


^G'G' 


{K) (K + G"; d) . (42) 


Substituting the formula of the intervalley dielectric ma¬ 
trix from Eq. (38), we obtain the following, remarkably 


simple analytical formula for the valley relaxation rate: 


Fv = 


riiEp 


(K + G) 


Vi{K + G;d) 


1 2 


ex 


(43) 


Equation (43) is the last one of the three key results of 
this work. 

Note that the screened result (43) can be ob¬ 


tained from the unscreened result (27) by substituting 


Vi {K -I- G; d) jf-K for V {K + G;d). On the one hand, 
it is remarkable that the matrix character of the inverse 


dielectric matrix in Eq. (38) does not appear explicitly 


in Eq. (43); instead, the effect of screening on the val¬ 


ley relaxation rate is described by a single scalar ck in 
Eq. (43). On the other hand, we emphasise that the 


role played by the quantity ck in the intervalley scatter¬ 
ing rate is analogous to the role played by the dielectric 
constant in the screening of a long-wave-length potential 
in a dielectric material. Accordingly, ejf can be called 
the effective intervalley dielectric constant. Importantly, 
the derivations of Eq. (27) and Eq. (43) rely on our spe¬ 
cific model for the electronic wave function (LCAO wave 
functions built from 2D Gaussian atomic wave functions) 
and the specific type of disorder (random, uncorrelated 
Goulomb impurities); therefore, it is possible that the 
notion of the effective intervalley dielectric constant is 
restricted to the present model only. 

In Fig.|^ the red curve shows the valley relaxation time 
for screened Coulomb impurities, according to Eq. (43), 
as a function of the effective Bohr radius UeB, for a 
graphene-impurity distance d = 0. The G sum of 


Eq. (43) was evaluated numerically. The valley relax¬ 


ation time is much longer in the screened case (red) than 
in the unscreened case (blue): apparently, screening is ef¬ 
fective in weakening the Coulomb potential of the impuri¬ 
ties even at the intervalley wave vector K, and therefore 
significantly prolongs the valley relaxation time. 
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Figure]^ also shows that for screened impurities, Tv is 
decreasing with increasing effective Bohr radius, whereas 
for unscreened impurities Tv shows an opposite trend. 
As the only difference between the corresponding results 
and (43) is the appearance of ek in the latter, the 


different trends are explained by the relatively fast decay 
of ek with increasing effective Bohr radius. 

At this point, the relative importance of the diagonal 
and off-diagonal matrix elements of the inverse interval¬ 
ley dielectric matrix for the screened result (red) of Fig.[^ 
is not known. Using the hypothesis that only the diago¬ 
nal matrix elements of e~^{K) are important, we calcu¬ 
late the valley relaxation rate (42) of a screened impurity 
with an ‘artificial’, diagonal inverse intervalley dielectric 
matrix, whose diagonal (off-diagonal) elements are given 
by e'^q{K) of Eq. (43) (are zero). The obtained val¬ 
ley relaxation rate is shown in Fig. |^as the green curve. 
The apparent qualitative difference between the screened 
result (red) and the diagonally screened result (green) re¬ 
veals that the off-diagonal matrix elements of the inverse 
intervalley dielectric matrix do play an important role in 
the screening of the charged impurities. 

Finally, Fig. shows the valley relaxation rate due 
to screened impurities as a function of the graphene- 
impurity distance, for three different values of the ef¬ 
fective Bohr radius. The behavior is very similar to the 
unscreened case, see Fig. [ijand the corresponding discus¬ 
sion in Sec. IIV Al 
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V. DISCUSSION AND CONCLUSIONS 

(1) In this work, we described the effects of screening 
using a model for the tt electrons of graphene, and we 
applied the linear-response framework of the RPA. It is 
a relevant, and, to our knowledge, open question how 
strongly the other (e.g., cr) bands influence the dielectric 
matrix at wave vector K. The substrate might also con¬ 
tribute to the screening of the short-wavelength Fourier 
components of the Coulomb scatterers. It is also an in¬ 
teresting future direction to describe how the effects be¬ 
yond linear-response behaviour, e.g^ bound-state forma¬ 
tion around the Coulomb impurityinfluence the val¬ 
ley relaxation time. In a related recent theory work, in¬ 
tervalley scattering due to a combined long-range-short- 
range scatterer was studied, using a method where the ef¬ 
fect of the long-range potential component was described 
nonperturbatively.^ 

(2) In order to simplify calculations and to reveal the 
role of the atomic structure of the electronic wave func¬ 
tion in the valley relaxation process, we have used 2D 
Gaussian wave functions [Eq. @] as building blocks of 
our model. This choice provided two advantages: (i) 
The 2D character of these atomic wave functions al¬ 
lows for a 2D description of screening effects, (ii) The 
Gaussian character allows one to derive simple analyti¬ 
cal results. We wish to emphasize that advantage (i) is 
rather substantial. Without (i), the description of screen- 


FIG. 6. Valley relaxation rate as a function of graphene- 
impurity distance (screened impurities). See caption of Fig.[^ 
for more details. 


ing would become much more involved: graphene has 
discrete translation invariance within its plane, but has 
no translational invariance in the out-of-plane direction. 
Hence, if 3D wave functions are used, then a ‘hybrid’ the¬ 
ory should be developed for screening, which would then 
incorporate a special dielectric linear response function 
E{q -\- G,—q ~ G',qz,q'z) [cf. Eq. (31)], where G and G' 
are 2D reciprocal lattice vectors, q is a wave vector from 
the 2D Brillouin zone, and q^ and q'^ are out-of-plane 
wave numbers. Developing such a hybrid theory would 
be a welcome development, but we do not attempt that 
in the present work. Advantage (ii) is not that substan¬ 
tial. In fact, all our results expressed with the form factor 
P{q) hold for any other 2D atomic wave function as well, 
as long as the latter is cylindrically symmetric. 

(3) Here, we described valley relaxation due to scat¬ 
tering off long-range Coulomb impurities. In general, the 
valley relaxation time is set by the interplay of a number 
of mechanisms (e.g., electron-phonon scattering, scatter¬ 
ing off short-range impurities, etc.). 

(4) Our calculation is restricted to zero temper¬ 
ature. At finite temperature, the valley relaxation 
rate is expected to change. One possible reason 

























for a temperature-dependent valley relaxation rate is 
temperature-dependent screening: the temperature- 
dependent electronic distribution /(e) appears in the for¬ 
mula for the polarizability matrix, see Eqs. (351 and (36). 
Another possible reason causing temperature-dependent 
Tv is the thermal population of phonons with large wave 
vectors, that are capable to scatter electrons between the 
valleys upon being absorbed. 

(5) Here we assumed spatial homogeneity of the elec¬ 
tronic distribution function /, and used the Boltzmann 
equation to describe the relaxation dynamics of a valley- 
polarized initial state. A complementary task is to de¬ 
scribe the valley dynamics in a spatially inhomogeous 
structure, where, e.g., valley-polarized electrons are in¬ 
jected to a nanostructure from a localised source.^! We 
expect that our Boltzmann-equation-based approach can 
be used as a starting point in that case, to derive macro¬ 
scopic transport equations describing valley diffusion, 
in analogy to the spin-diffusion equations developed in 
spintronics.l^ 

(6) The robustness of the valley index against scat¬ 
tering processes can be characterised by the valley-flip 
length or intervalley scattering length Lp. this is the typ¬ 
ical distance an electron can travel without having its 
valley index flipped. In our case of charged impurities, 
intravalley scattering happens more often than interval¬ 
ley scattering, hence the motion of the electron between 
two valley flips is diffusive. Our results for the valley re¬ 
laxation time Tv allow us to estimate the dependence of 
the valley-flip length as a function of system parameters. 
The diffusion coefficient is estimated as D = Ttru|/2, 
where Ttr is the transport lifetime. Then the valley-flip 
length is Li = y/TJri, where t; = 2tv is the interval¬ 
ley scattering time, see Eq. 0- As Ttr oc e-pjni (see 
Ref. [5D| and t „ oc l/me-p [see Eqs. (27) and (43)], the 
valley-flip length is independent of the Fermi energy ep 
and inversely proportional to the impurity sheet density 
TT-i. In Fig. we show the valley-flip length Li as a 
function of the impurity sheet density and the graphene- 
impurity distance d, with the the effective Bohr radius 
set to OeB = 0.2acc; to obtain this result, Eq. (3.23) of 
Ref. [Sni was used for Ttr and our result (431 was used for 
Tv. Note that besides impurities, the rough edges of a real 
graphene sample can also induce intervalley scattering.!^ 
The relative importance of edge-induced and Coulomb- 
impurity-induced intervalley scattering can be judged by 
comparing the sample size and the valley-flip length eval¬ 
uated in Fig. For example, for a sample size of 10 mi¬ 
crons and graphene-impurity distance of d = accj Fig-0 
suggests that Coulomb scattering gains importance over 
edge scattering if rii > 3 x 10^^ cm“^. 

(7) In a recent measurement of the valley Hall effect 
in graphene,the length scale characterising the spa¬ 
tial decay of the non-local resistance was found to be 
« 1.0 /im. This length scale can be identifiecP^ as the 
valley-flip length Li we defined above. Charged impu¬ 
rities are thought to be present in the measured device 
(see section 6 of the Supplementary Material of Ref. nil). 
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FIG. 7. Valley-flip length as a function of the impurity sheet 
density and graphene-impurity distance. The effective Bohr 
radius is set to OeB = 0.2acc. The unit of horizontal axis is 
the carbon-carbon distance acc. 



hence it is motivated to relate this measured length scale 
to our theoretical results. Assuming that the charged im¬ 
purity sheet density in the sample is in the range shown 
in our Fig. 7, the following interpretations can be sug¬ 
gested: the measured length scale Li « 1 ^m is set by (i) 
charged impurities that are very close to the graphene 
plane {d < 0.5acc), or (ii) sources other than charged 
impurities, e.g., edges or short-range impurities. 

In conclusion, we have presented a model for valley 
relaxation due to randomly positioned charged impuri¬ 
ties in graphene. We described the dependence of the 
valley relaxation rate of an ensemble of valley-polarized 
electrons on the model parameters (Fermi energy, im¬ 
purity sheet density, graphene-impurity distance, spatial 
extension of the atomic Pz wave functions). The static 
screening of the charged impurities was described by, as 
required for crystalline materials, the dielectric matrix, 
which we evaluated in the RPA. Our results highlight 
that a quantitatively accurate description of valley re¬ 
laxation is more challenging than that of the electrical 
conductivity: the former requires that screening due to 
electron-electron interaction is described in terms of the 
dielectric matrix, and that the spatial variation of the 
electronic wave functions on the atomic length scale is 
precisely known. 
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Appendix A: Conventions 

In this Appendix, the reference frame is specified and 
the vectors characterising the direct and reciprocal lat¬ 
tices are defined. The direct lattice is shown in Fig. 
Atoms of the A and B sublattices are depicted as black 
points and circles, respectively. The shaded rhombus 
shows the unit cell of the direct lattice. The primitive 
vectors of the direct lattice are 


ai = 


acc f -a/3 


02 = 


QCC /Vs 

2 I 3 


(Al) 


and the vector connecting the A and B sites within a unit 
cell is 


T = acc 


The primitive vectors of the reciprocal lattice are 


(A2) 


bi = 


27T 


3acc 

The Dirac points are 
47r 


3,1 b2 = 


27r 

3acc 


(-V3, l) . (A3) 


K = 


^ 3v^acc 


Appendix B: Proof of Eqs. (8a| and (8bl 


First, we prove Eq. (8aI, which corresponds to the in¬ 


travalley transition probability. The formula of the ma¬ 
trix element Vc,k+k,c,k+k’ can be approximated by keep¬ 
ing only those terms that are proportional to V{Ak), 
and dropping all other terms involving V (Ak + G) with 
G / 0 terms. This is a reasonable approximation as our 
small-Fermi-energy condition implies V{Ak) ~ 1/Ak ^ 
1/G ~ V{Ak + G). With this simplification, we find 


v; 


C,Jir+K,C,JC + K,' 


-Sc,K+n,c,K+H.' (Ak) P (Ak) V* (Ak) , 

(Bl) 

A similar formula holds for valley K' with 
Sc.k'+k.,c,k'+k' ■ The transition probabilities are 
equal in the valleys if the structure factor has the 
property. 


\Sc^K+K,,c,K+n' (^^)| — \Sc.K'+ k.,c,K '(^^)l 




which can be proven by substituting Eqs. (17a I and (17b) 
to the definition of structure factor and approximat- 
1 due to Ak —^ 0. This argument can be 


mg e 


— iA« 


reused for the case of screened impurities, assuming that 
long-wave-length screening is appropriately descr ibed by 
the frequently used jellium-type RPA description.l^nMESl 


To prove Eq. (8b I, which corresponds to the intervalley 


transition rates, we rewrite Eq. (8b) as 


Wk+k,K'+k' — Wk+k.',K'+k — Wk'+k.,K+k.', (B3) 

where the second equality expresses detailed balance, 
which arises here as the transition rates are evaluated 
from Fermi’s Golden Rule and hence are invariant for 
the exchange of the initial and the final states. Thus, 
to confirm the first equality, it is enough to investi¬ 
gate the dependence of the intervalley transition rate 
Wk+k..k'+k.' on the angles (p and ip' of k and k', re¬ 
spectively. The transition rate depends on ip and ip' only 
through jc'+k' (-^ + which is invariant un¬ 

der the exchange of ip and (/p'polar angle of k and n' 
as seen from Eq. (23). This proves the first equality 


in Eq. (B3), and hence (8b), for both unscreened and 


screened impurities. 


Appendix C: Relation of the dielectric matrix and 
the polarizability matrix 


In this Appendix, we establish the relation (32) be¬ 


tween the dielectric matrix and the polarizability matrix. 

First, recall that the real-space dielectric function 
e(r,r') describes the relation between the external and 
total potentials [cf. Eq. (29)], 


V'ext(r) = J (fr'e{r,r')Vtotir'), (Cl) 

whereas the real-space polarizability function n(r, r') de¬ 
scribes the relation between the induced electron density 
and the total potential. 


n-ind (r) = J d^r'n (r,r') Ftot (r'). 


(C2) 


For our purposes, r, r' are 2D position vectors, since the 
model we use to describe graphene’s electrons is 2D. 

Next, we apply Fourier transformation on the defini¬ 
tions and ( |C2[ ), exploiting their invariance with re¬ 
spect to lattice translations: e.g., H (r-|-i?, r'-|-i?) = 
n (r, r') for any lattice vector R. Recall that we use the 
following definition for Fourier transformation 


f{k)=J d^rf{r)< 


— ikr 


(C3) 


where f{r) is an arbitrary 2D position-dependent func¬ 
tion. We evaluate the Fourier transforms of Eqs. (Cl) 


and (C2) at an arbitrary wave vector k = q + G, where q 


is within the Brillouin zone and G is a reciprocal lattice 
vector; we obtain 

Eext {q + G)=J2 ecG' (g) ^tot {q + G'), (C4) 


G' 


“^ind 


iq + G)=J2^GG'iq)Vtot{q + G'), (C5) 


G' 
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FIG. 8 . Direct and reciprocal lattices, (a) Direct lattice of graphene, showing the unit cell (shaded rhombus), the A (points) 
and B (circles) sublattices, the primitive lattice vectors ai and 02 , the vector r connecting the A and B sites within a unit 
cell, and the reference frame, (b) Reciprocal lattice of graphene (black points), showing the Brillouin zone (shaded hexagon), 
the primitive reciprocal lattice vectors 61 and b 2 , the special points K and K' of the Brillouin zone, the shifted reciprocal 
lattices consisting of points K + G and K' + G (blue upward triangles and red downward triangles, respectively). The dashed 
circle shows the reciprocal-space extension of the Gaussian atomic electron density, i.e., the set of momentum vectors where 
the Fourier-transformed atomic electron density is P{k) = 1/2, for an effective Bohr radius of OeB = 0.2acc. 


Here, we used the fact that the invariance of the real- 
space response function with respect to lattice transla¬ 
tions implies, e.g., n(q + G,q' + G') oc (5q,-q' and in¬ 
troduced the dielectric matrix e(< 7 ) and the polarizabil¬ 
ity matrix n(q) as cgg' (?) = e (9 + G, —q — G') and 
ncG' (g) = n (q -b G, -q - G'), respectively. 

Having these definitions at hand, the dielectric ma¬ 
trix and the polarizability matrix can be connected using 
Coulomb’s law. Coulomb’s law implies that the induced 
potential, i.e., the potential created by the induced elec¬ 
tron density, reads 

l^ind {r) = J d^r' Vb (r - r') mnd {r'), (C6) 

where Vc {r — r') = ep/jr — r'| is the 2D Coulomb po¬ 
tential. Applying Fourier transform with respect to r 
yields 


Mnd (q + G) — Vc (q + G) riind (q + G) , 


(C7) 


where Vc {q + G) = 2Trel/\q+G\ is the 2D Fo urier trans¬ 
form of the Coulomb potential. Inserting Eq. (C71 to the 
relation 


Hot (q + G) — I4xt (q + G) + Hnd (q + G) 


(C8) 


and eliminating nind via Eq. (C5), one finds 


Hxt (9 + G)=i: [<^GG' - He (q + G) Hgg' ( 9 )] 

G' 

xHot(q + G'), (C9) 


proving Eq. (32). 


Appendix D: Estimation of the effective Bohr radius 


In the main text, the atomic Pz orbital is represented 
by the 2D Gaussian-type wave function </>('*')) defined in 
Eq. (131. The 2D spatial extension of this wave function 


in the xy plane is characterized by the effective Bohr 
radius UeB- Here, we estimate the 2D spatial extension of 
a three-dimensional 2pz orbital of a free-standing carbon 
atom, to provide an estimate for OeB- 

To this end, we invoke those results of Ref. [53] that 
correspond to a free-standing carbon atom. In Ref. 1631 
the three-dimensional single-electron orbitals are approx¬ 
imated by Slater-type orbitals; in particular, the 2pz or¬ 
bital of a carbon atom is approximated by 


(t>pA'r,z) = 


■\/327r(ao7^effF 


_ vH+ii 

e 2ao/^eff 


(Dl) 


where r = {x,y), r = ap = 0.53A is the 

Bohr radius, and Zes = 3.14 is the effective charge of the 
nucleus. (The latter is denoted in Ref. [531 as Z — a.) 

We characterize the 2D spatial extension of (j>p^ in the 
xy plane via 


(bpjr IbpJ =——ap/Zeff « O.SOA (D2) 
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The same quantity for our 2D Gaussian-type wave func- Expressing this result in units of the carbon-carbon dis- 
tion is tance acc, we find OgB ~ 0.28acc- 

(D3) 

From the requirement {(pp^ \ r \4>p^) = {4>\ t \4‘)^ we obtain 
the estimate 
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aeB-y 


ao/^eff ~ 0-40 A. 


(D4) 


^ A. Rycerz, J. Tworzydlo, and C. W. J. Beenakker, Nat. 
Phys. 3, 172 (2007). 

^ D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 
(2007). 

® W. Yao, D. Xiao, and Q. Nin, Phys. Rev. B 77, 235406 
(2008). 

^ D. Gunlycke and C. T. White, Phys. Rev. Lett. 106, 
136806 (2011). 

® W.-K. Tse, A. Saxena, D. L. Smith, and N. A. Sinitsyn, 
Phys. Rev. Lett. 113, 046602 (2014). 

® W.-Y. Shan, J. Zhou, and D. Xiao, Phys. Rev. B 91, 
035402 (2015). 

^ L. E. Golub and S. A. Tarasenko, Phys. Rev. B 90, 201402 
(2014). 

® T. O. Wehling, A. Huber, A. 1. Lichtenstein, and M. 1. 

Katsnelson, Phys. Rev. B 91, 041404 (2015). 

® Z. Zhu, A. Gollaudin, B. Fauque, W. Kang, and K. Behnia, 
Nat. Phys. 8, 89 (2012). 

“ T. Gao, G. Wang, W. Han, H. Ye, G. Zhu, J. Shi, Q. Niu, 
P. Tan, E. Wang, B. Liu, et al., Nat. Gommun. 3, 887 
( 2012 ). 

K. F. Mak, K. He, J. Shan, and T. F. Heinz, Nat. Nano. 
7, 494 (2012). 

H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Nat. Nano. 
7, 490 (2012). 

J. Isberg, M. Gabrysch, J. Hammersberg, S. Majdi, K. K. 
Kovi, and D. J. Twitchen, Nat. Mater. 12, 760 (2013). 

E. A. Laird, F. Pei, and L. P. Kouwenhoven, Nat. Nan¬ 
otech. 8, 565 (2013). 

R. V. Gorbachev, J. C. W. Song, G. L. Yu, A. V. Kretinin, 

F. Withers, Y. Gao, A. Mishchenko, 1. V. Grigorieva, K. S. 
Novoselov, L. S. Levitov, et al., Science 346, 448 (2014). 

K. F. Mak, K. L. McGill, J. Park, and P. L. McEuen, 
Science 344, 1489 (2014). 

P. Recher, B. Trauzettel, A. Rycerz, Y. M. Blanter, 

C. W. J. Beenakker, and A. F. Morpurgo, Phys. Rev. B 
76, 235404 (2007). 

A. Palyi and G. Burkard, Phys. Rev. Lett. 106, 086801 

( 2011 ). 

G. Y. Wu, N.-Y. Lue, and L. Chang, Phys. Rev. B 84, 
195463 (2011). 

G. Y. Wu and N.-Y. Lue, Phys. Rev. B 86, 045456 (2012). 
G. Y. Wu, N.-Y. Lue, and Y.-C. Chen, Phys. Rev. B 88, 
125422 (2013). 

D. Culcer, A. L. Saraiva, B. Koiller, X. Hu, and 

S. Das Sarma, Phys. Rev. Lett. 108, 126804 (2012). 

G. Szechenyi and A. Palyi, Phys. Rev. B 89, 115409 (2014). 
A. Kormanyos, V. Zolyomi, N. D. Drummond, and 
G. Burkard, Phys. Rev. X 4, 011034 (2014). 


G. -B. Liu, H. Pang, Y. Yao, and W. Yao, New Journal of 
Physics 16, 105011 (2014). 

®® N. Rohling and G. Burkard, New Journal of Physics 14, 
083008 (2012). 

N. Rohling, M. Russ, and G. Burkard, Phys. Rev. Lett. 
113, 176801 (2014). 

A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007). 
2® A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 
(2009). 

T. Ando, Journal of the Physical Society of Japan 75, 
074716 (2006). 

E. H. Hwang, S. Adam, and S. D. Sarma, Phys. Rev. Lett. 
98, 186806 (2007). 

Y. Song and H. Dery, Phys. Rev. Lett. Ill, 026601 (2013). 

H. Ochoa, F. Finocchiaro, F. Guinea, and V. 1. Fal’ko, 
Phys. Rev. B 90, 235429 (2014). 

C. Mai, Y. G. Semenov, A. Barrette, Y. Yu, Z. Jin, L. Gao, 
K. W. Kim, and K. Gundogdu, Phys. Rev. B 90, 041414 
(2014). 

A. Pachoud, A. Ferreira, B. Ozyilmaz, and A. H. Cas¬ 
tro Neto, Phys. Rev. B 90, 035444 (2014). 

®® L. Braginsky and M. Entin, arXiv:1412.7810 (unpub¬ 
lished) . 

G. Csiszar and A. Palyi, Phys. Rev. B 90, 245413 (2014). 
C. Tahan and R. Joynt, Phys. Rev. B 89, 075302 (2014). 

®® C. H. Yang, A. Rossi, R. Ruskov, N. S. Lai, F. A. Mo- 
hiyaddin, S. Lee, C. Tahan, G. Klimeck, A. Morello, and 
A. S. Dzurak, Nat. Gommun. 4 (2013). 

S. L. Adler, Phys. Rev. 126, 413 (1962). 

N. Wiser, Phys. Rev. 129, 62 (1963). 

M. van Schilfgaarde and M. 1. Katsnelson, Phys. Rev. B 
83, 081409 (2011). 

T. Tudorovskiy and S. A. Mikhailov, Phys. Rev. B 82, 
073411 (2010). 

H. Suzuura and T. Ando, Phys. Rev. Lett. 89, 266603 

( 2002 ). 

E. McCann, K. Kechedzhi, V. 1. Fal’ko, H. Suzuura, 
T. Ando, and B. L. Altshuler, Phys. Rev. Lett. 97, 146805 
(2006). 

‘‘® S. V. Morozov, K. S. Novoselov, M. 1. Katsnelson, 

F. Schedin, L. A. Ponomarenko, D. Jiang, and A. K. Geim, 
Phys. Rev. Lett. 97, 016801 (2006). 

X. Wu, X. Li, Z. Song, C. Berger, and W. A. de Heer, 
Phys. Rev. Lett. 98, 136801 (2007). 

F. V. Tikhonenko, D. W. Horsell, R. V. Gorbachev, and 
A. K. Savchenko, Phys. Rev. Lett. 100, 056802 (2008). 

‘^® L. Malard, M. Pimenta, G. Dresselhaus, and M. Dressel- 
haus. Physics Reports 473, 51 (2009). 



15 


D. Graf, F. Molitor, K. Ensslin, C. Stampfer, A. Jungen, 

C. Hierold, and L. Wirtz, Nano Letters 7, 238 (2007). 

P. R. Wallace, Phys. Rev. 71, 622 (1947). 

S. Das Sarnia, S. Adam, E. H. Hwang, and E. Rossi, Rev. 
Mod. Phys. 83, 407 (2011). 

E. Fradkin, Phys. Rev. B 33, 3257 (1986). 

E. Fradkin, Phys. Rev. B 33, 3263 (1986). 

E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 
(2007). 

B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Jour¬ 
nal of Physics 8, 318 (2006). 

J. Lindhard and A. Winther, Mat. Fys. Medd. Dan. Vid. 


Selsk. 34, 4 (1964). 

V. M. Pereira, J. Nilsson, and A. H. Gastro Neto, Phys. 
Rev. Lett. 99, 166802 (2007). 

A. V. Shytov, M. 1. Katsnelson, and L. S. Levitov, Phys. 
Rev. Lett. 99, 236801 (2007). 

D. S. Novikov, Phys. Rev. B 76, 245435 (2007). 

T. Valet and A. Fert, Phys. Rev. B 48, 7099 (1993). 

D. A. Abanin, A. V. Shytov, L. S. Levitov, and B. 1. 
Halperin, Phys. Rev. B 79, 035304 (2009). 

E. Glementi and D. L. Raimondi, The Journal of Ghemical 
Physics 38, 2686 (1963). 



